9. Γεωεπεξεργασία διανυσματικών δεδομένων
Περιεχόμενα
9. Γεωεπεξεργασία διανυσματικών δεδομένων¶
Ειδική ενότητα για εκτέλεση στο Google Colab¶
# έλεγχος αν το notebook τρέχει στο google colab
try:
import google.colab
IN_COLAB = True
except:
IN_COLAB = False
# αν το notebook τρέχει στο colab, mount το Google Drive και αλλαγή στο directory που έχει γίνει clone το github repository.
# εγκατάσταση απαραίτητων βιβλιοθηκών
if IN_COLAB:
from google.colab import drive
drive.mount('/content/drive')
%cd /content/drive/MyDrive/Colab\ Notebooks/programming/notebooks
!pip install rasterio geopandas folium matplotlib mapclassify
Ο χώρος και οι γεωμετρικές δομές του μπορούν να αναπαρασταθούν μέσω διανυσματικών δεδομένων (Vector). Οι τρείς βασικές δομές δεδομένων είναι:
Τα σημεία
Οι γραμμές
Τα πολύγωνα

Επιπλέον αυτά τα γεωμετρικά δεδομένα συνοδεύονται και απο περιγραφικά δεδομένα που αφορούν τις ιδιότητες ή τα χαρακτηριστικά αυτών των δεδομένων.

Στα Σ.Γ.Π. ο πιο συνηθισμένος τύπος αρχείων αποθήκευσης αυτών των δεδομένων είναι το shapefile. Πλέον έχουν αναπτυχθεί και άλλοι τύποι όπως geojson, geopackage και χωρικές βάσεις δεδομένων (geodatabases) όπως η Postgresql/Postgis.
Ο προσδιορισμός της γεωγραφικής θέσης στην υδρόγειο γίνεται μέσω ενός ζεύγους γεωγραφικών συντεταγμένων. Κάθε σημείο στον χώρο προσδιορίζεται γεωγραφικά από το γεωγραφικό μήκος (λ) και το γεωγραφικό πλάτος (φ).

Για να αποδοθεί η τρισδιάστατη υδρόγειος σφαίρα σε ένα δυσδιάστατο σύστημα αναφορά χρησιμοποιείται ένα προβολικό σύστημα.

Κάθε προβολικό σύστημα της σφαίρας στο επίπεδο εισάγει μια σειρά παραμορφώσεων που αφορά το σχήμα των γεωμετρικών δομών, την κλίμακα, την έκταση και τις αποστάσεις. Ανάλογα το προβολικό σύστημα κάποιες από τις παραπάνω παραμορφώσεις εμφανίζονται σε μεγάλο βαθμό και άλλες όχι. Οπότε ανάλογα το είδος της έρευνας ο ερευνητής οφείλει να γνωρίζει τι είδους παραμορφώσεις εισάγει η κάθε προβολή και ανάλογα να επιλέγει την προβολή με τα λιγότερα σφάλματα.
Python βιβλιοθήκες για διανυσματικά δεδομένα¶
Για την ανάγνωση, εγγραφή, επεξεργασία διανυσματικών δεδομένων στην Python έχουν καθιερωθεί μια σειρά βιβλιοθηκών. Η αρχαιότερη και βασική βιβλιοθήκη είναι η GDAL/OGR. Επειδή η βιβλιοθήκη δεν είναι ιδιαίτερα συμβατή σε σχέση με τον τρόπο συγγραφής της Python και είναι επιρρεπής σε σφάλματα έχουν αναπτυχθεί πιο σύγχρονες βιβλιοθήκες όπως η βιβλιοθήκη Fiona που είναι ιδιαίτερα χρήσιμη για την ανάγνωση/εγγραφή διανυσματικών δεδομένων και η βιβλιοθήκη Shapely η οποία χρησιμοποιείται για επεξεργασία και ανάλυση δεδομένων.
Η βιβλιοθήκη Shapely¶
shapely from WKT¶
Μπορούμε να δημιουργήσουμε αντικείμενα shapely που αναπαριστούν σημεία ή γραμμές ή πολύγωνα μέσω WKT. Η γλώσσα WKT είναι μια ειδική διάλεκτος για την περιγραφή διανυσματικών αντικειμένων. Εισάγουμε τις απαραίτητες υπο-βιβλιοθήκες (geometry και wkt) από την βιβλιοθήκη shapely.
import shapely.geometry
import shapely.wkt
Πολύγωνα¶
Καλούμε την μέθοδο shapely.wkt.loads() για να δημιουργήσουμε shapely objects από wkt
pol1 = shapely.wkt.loads("POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0,0 0))")
pol1
Η εκτύπωση του shapely αντικειμένου επιστρέφει μία περιγραφή σε μορφή WKT.
print(pol1)
POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0))
Εκτύπωση του τύπου του αντικειμένου pol1
type(pol1)
shapely.geometry.polygon.Polygon
pol2 = shapely.wkt.loads("POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))")
print(pol2)
POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))
pol2
Δημιουργία MultiPolygon shapely object από αντίστοιχη συλλογή πολυγώνων MULTIPOLYGON
pol3 = shapely.wkt.loads("""
MULTIPOLYGON
(((40 40, 20 45, 45 30, 40 40)),
((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))
""")
type(pol3)
shapely.geometry.multipolygon.MultiPolygon
pol3
type(pol3)
shapely.geometry.multipolygon.MultiPolygon
Σημεία¶
pnt1 = shapely.wkt.loads("""
POINT (30 10)
""")
pnt1
Συλλογή σημείων
pnt2 = shapely.wkt.loads("""
MULTIPOINT(0 0,1 1)
""")
pnt2
Γραμμές¶
line1 =shapely.wkt.loads("""
LINESTRING(1.5 2.45,3.21 4)
""")
line1
Συλλογή γραμμών
line2 =shapely.wkt.loads("""
MULTILINESTRING((0 0,-1 -2,-3 -4, -3 -8),(2 3,3 4,6 7))
""")
line2
Shapely μέσω συναρτήσεων¶
from shapely.geometry import Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon
Κάθε συνάρτηση λαμβάνει διαφορετικά ορίσματα ανάλογα με την συνάρτηση.
Στην συνέχεια ακολουθεί η δημιουργία ενός αντικειμένου shapely γεωμετρίας σημείου.
pnt1 = shapely.geometry.Point((2, 0.5))
pnt1
Για την δημιουργία συλλογής σημείων multipoint δίνουμε σαν όρισμα στην συνάρτηση shapely.geometry.MultiPoint μία λίστα από πλειάδες (tuples). Η κάθε πλειάδα (x,y)αντιστοιχεί στις συντεταγμένες ενός σημείου.
coords = [(2, 0.5), (1, 1), (-1, 0), (1, 0)]
pnt2 = shapely.geometry.MultiPoint(coords)
pnt2
Για την δημιουργία γραμμών χρησιμοποιείται πάλι μία λίστα από tuples που περιγράφουν τις κορυφές της γραμμής. Στο παρακάτω παράδειγμα χρησιμοποιούμαι την προηγούμενη πλειάδα αλλά πλέον χρησιμοποιούμε την μέθοδο shapely.geometry.LineString για την δημιουργία γραμμής και όχι την μέθοδο shapely.geometry.MultiPoint που δημιουργεί συλλογές σημείων.
line1 = shapely.geometry.LineString(coords)
line1
Κατασκεύη γραμμής από Point Objects:
p1 = shapely.wkt.loads("POINT (0 0)")
p2 = shapely.wkt.loads("POINT (1 1)")
p3 = shapely.wkt.loads("POINT (2 -1)")
p4 = shapely.wkt.loads("POINT (2.5 2)")
p5 = shapely.wkt.loads("POINT (1 -1)")
shapely.geometry.LineString([p1, p2, p3, p4, p5])
Αντίστοιχα μπορούμε να δημιουργήσουμε συλλογές γραμμών.
Δημιουργούμε ξεχωριστά αντικείμενα γραμμών shapely και στην συνέχεια καλούμε την συνάρτηση shapely.geometry.MultiLineString όπου ορίζουμε σαν όρισμα μια λίστα με τις μεμονωμένες γραμμές.
l1 = shapely.geometry.LineString([(2, 0.5), (1, 1), (-1, 0), (1, -1)])
l2 = shapely.geometry.LineString([(-2, 1), (2, -0.5), (3, 1)])
line2 = shapely.geometry.MultiLineString([l1, l2])
line2
Αντίστοιχα δημιουργούμε πολυγωνικά αντικείμενα μέσω μια λίστας με πλειάδες σημείων που χρησιμοποιείται σαν όρισμα στην συνάρτηση shapely.geometry.Polygon
coords = [(0, 0), (0, -1), (7.5, -1), [7.5, 0], (0, 0)]
shapely.geometry.Polygon(coords)
Αν θέλουμε μπορούμε να περάσουμε μία δεύτερη λίστα η οποία περιλαμβανει επιμέρους λίστες με πλειάδες σημείων τα οποία περιγράφουν τρύπες μέσα στο πολύγωνο.
coords_exterior = [(0, 0), (0, -1), (7.5, -1), (7.5, 0), (0, 0)]
coords_interiors = [[(0.4, -0.3), (5, -0.3), (5, -0.7), (0.4, -0.7), (0.4, -0.7)]]
shapely.geometry.Polygon(coords_exterior, coords_interiors)
Με την συνάρτηση shapely.geometry.MultiPolygon δημιουργούμε αντίστοιχα συλλογές πολυγώνων από μεμονωμένα αντικείμενα shapely.geometry.polygon.Polygon
pol2
multipolygon1 = shapely.geometry.MultiPolygon([pol1, pol2])
multipolygon1
print(multipolygon1)
MULTIPOLYGON (((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0)), ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1)))
Σύμφωνα με τις προδιαγραφές simple features το παραπάνω object δεν είναι έγκυρο γιατί ένα πολύγωνο τέμνει ένα άλλο σε άπειρο αριθμό σημείων. Μπορούμε να ελέγξουμε την εγκυρότητα ενός αντικειμένου με την κλήση τις ιδιότητας is_valid. Επίσης κατά την οπτικοποίηση στο προηγούμενο βήμα του αντικειμένου multipolygon1 αυτό εμφανίζεται με κόκκινο (αντί για πράσινο). Ένδειξη ότι δεν ειναι valid.
multipolygon1.is_valid
False
Ταυτόχρονα μπορούμε να φτιάξουμε σύνθετες συλλογές από επιμέρους αντικείμενα shapely.
geo_collection = shapely.geometry.GeometryCollection([multipolygon1,line2,pnt1])
geo_collection
Shapely μέσω shape¶
Μπορούμε να δημιουργήσουμε αντικείμενα shapely μέσω της συνάρτησης shapely.geometry.shape η οποία δέχεται σαν όρισμα ένα λεξικό μορφής GEOJSON το οποίο πρέπει να έχει τις εξής δύο ιδιότητες.
Την ιδιότητα
"type"που περιγράφει τον τύπο γεωμετρίαςΤην ιδιότητα
"coordinates"που περιγράφει τις γεωμετρίες και τις συντεταγμένες τους σαν λίστες ή πλειάδες.
d = {"type": "Point", "coordinates": (0, 1)}
shapely.geometry.shape(d)
Δημιουργία αντικειμένου MultiPolygon μέσω GEOJSON λεξικού:
d = {
"type": "MultiPolygon",
"coordinates": [
[
[[40, 40], [20, 45], [45, 30], [40, 40]]
],
[
[[20, 35], [10, 30], [10, 10], [30, 5], [45, 20],
[20, 35]],
[[30, 20], [20, 15], [20, 25], [30, 20]]
]
]
}
pol3 = shapely.geometry.shape(d)
pol3
Γεωμετρικός τύπος δεδομένων¶
H ιδιότητα .geom_type ενός αντικειμένου shapely περιγράφει τον γεωμετρικό τύπο της:
pol1.geom_type
'Polygon'
line1.geom_type
'LineString'
geo_collection.geom_type
'GeometryCollection'
Συντεταγμένες¶
Για να ανακτήσουμε τις συντεταγμένες ενός shapely αντικειμένου καλούμε με διαφορετικό τρόπο τις ανάλογες συναρτήσεις ανάλογα την πολυπλοκότητα του κάθε τύπου.
Για ένα Point object καλούμε άμεσα την ιδιότητα coords η οποία επιστρέφει ένα shapely.coords.CoordinateSequence object.
pnt1.coords
<shapely.coords.CoordinateSequence at 0x7f5d340f3970>
Μπορούμε να λάβουμε ως λίστα τις πλειάδες συντεταγμένων που το απαρτίζουν
list(pnt1.coords)
[(2.0, 0.5)]
Αντίστοιχα για γραμμή
list(line1.coords)
[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, 0.0)]
Για να ελέγξουμε σε ένα αντικείμενο συλλογών γεωμετρίας (MultiPoint, MultiPolygon κτλ) πόσες επιμέρους γεωμετρίες περιλαμβάνει:
len(line2.geoms)
2
line2
Για να λάβουμε το πρώτο αντικείμενο γεωμετρίας:
line2.geoms[0]
type(line2.geoms[0])
shapely.geometry.linestring.LineString
line2.geoms[0].geom_type
'LineString'
Και λαμβάνουμε τις συντεταγμένες της πρώτης γεωμετρίας:
list(line2.geoms[0].coords)
[(2.0, 0.5), (1.0, 1.0), (-1.0, 0.0), (1.0, -1.0)]
ή της δεύτερης
list(line2.geoms[1].coords)
[(-2.0, 1.0), (2.0, -0.5), (3.0, 1.0)]
Για τα πολύγωνα ακολουθούμε διαφορετική προσέγγιση. Ένα πολύγωνο αποτελεί από το εξωτερικό περίγραμμα (exterior) ή και ένα ή περισσότερα περιγράμματα εσωτερικών οπών (interiors). Κατά συνέπεια έχουμε συντεταγμένες που περιγράφουν το κάθε περίγραμμα.
Το εξωτερικό περίγραμμα του pol1 αντικειμένου:
pol1.exterior
Και οι συντεταγμένες του:
list(pol1.exterior.coords)
[(0.0, 0.0), (0.0, -1.0), (7.5, -1.0), (7.5, 0.0), (0.0, 0.0)]
pol1
Το pol1 όμως δεν έχει εσωτερικές τρύπες γι αυτό και το παρακάτω επιστρέφει μηδέν.
len(pol1.interiors)
0
Έστω το παρακάτω MultiPolygon object που δημιουργήσαμε σε προηγούμενο στάδιο.
pol3
Περιλαμβάνει δύο ξεχωριστά γεωμετρικά αντικείμενα:
len(pol3.geoms)
2
Ας πάρουμε το πρώτο:
pol3.geoms[0]
Δεν περιλαμβάνει καμία τρύπα στο εσωτερικό του. Ας πάρουμε τις συντεταγμένες από το περίγραμμά του (exterior)
pol3.geoms[0].exterior.coords
<shapely.coords.CoordinateSequence at 0x7f5d34139a80>
Και ας τις επιστρέψουμε σαν λίστα πλειάδων από ζεύγη της μορφής (x,y)
list(pol3.geoms[0].exterior.coords)
[(40.0, 40.0), (20.0, 45.0), (45.0, 30.0), (40.0, 40.0)]
Ας δοκιμάσουμε το δεύτερο αντικείμενο γεωμετρίας. Ας το δούμε:
pol3.geoms[1]
Λαμβάνουμε τις συντεταγμένες του εξωτερικού περιγράμματος:
list(pol3.geoms[1].exterior.coords)
[(20.0, 35.0),
(10.0, 30.0),
(10.0, 10.0),
(30.0, 5.0),
(45.0, 20.0),
(20.0, 35.0)]
Ας δούμε πόσες τρύπες έχει στο εσωτερικό του:
len(pol3.geoms[1].interiors)
1
Παίρνουμε το περίγραμμα της τρύπας
pol3.geoms[1].interiors[0]
Και τις συντεταγμένες της
list(pol3.geoms[1].interiors[0].coords)
[(30.0, 20.0), (20.0, 15.0), (20.0, 25.0), (30.0, 20.0)]
Ιδιότητες αντικειμένων¶
Υπολογισμός ορίων (bounds)¶
geo_collection
geo_collection.bounds
(-2.0, -1.0, 7.5, 1.0)
shapely.geometry.box(*geo_collection.bounds)
line1
line1.bounds
(-1.0, 0.0, 2.0, 1.0)
pnt1
pnt1.bounds
(2.0, 0.5, 2.0, 0.5)
list(pnt1.coords)
[(2.0, 0.5)]
Υπολογισμός μήκους γραμμής¶
line1
line1.length
5.354101966249685
line2.geoms[0]
line2.geoms[0].length
5.5901699437494745
Νέες γεωμετρίες¶
Κεντροειδές πολυγώνου (centroid)¶
pol2
pol2.centroid
shapely.geometry.GeometryCollection([pol2, pol2.centroid])
Περιμετρική ζώνη (buffer)¶
pnt1.buffer(5)
shapely.geometry.GeometryCollection([pnt1,pnt1.buffer(5)])
pol1.buffer(5)
shapely.geometry.GeometryCollection([pol1,pol1.buffer(5)])
Σχέσεις μεταξύ αντικειμένων¶
shapely.geometry.GeometryCollection([pol1,pol3])
pol3
pol1.intersects(pol3)
False
shapely.geometry.GeometryCollection([pol1,pol2])
pol1.intersects(pol2)
True
Γεωμετρικές πράξεις¶
x = shapely.geometry.Point((0, 0)).buffer(1)
y = shapely.geometry.Point((1, 0)).buffer(1)
shapely.geometry.GeometryCollection([x, y])
x.intersection(y)
x.difference(y)
x.union(y)
x.union(y)
Υπολογισμός απόστασης ανάμεσα σε δύο αντικείμενα
shapely.geometry.GeometryCollection([pol1, pol3])
pol1.distance(pol3)
10.307764064044152
pol3.distance(pol1)
10.307764064044152
Μετασχηματισμός προβολικού συστήματος:
import pyproj
from shapely.geometry import Point
from shapely.ops import transform
wgs84_pt = Point(39.35858398397631, 22.932070043920394)
wgs84 = pyproj.CRS('EPSG:4326')
greek_grid = pyproj.CRS('EPSG:2100')
project = pyproj.Transformer.from_crs(wgs84, greek_grid, always_xy=True).transform
projected_point = transform(project, wgs84_pt)
print(projected_point)
POINT (2087909.8290560723 2620022.621513317)
list(wgs84_pt.coords)[0]
(39.35858398397631, 22.932070043920394)
Προβολή δεδομένων σε διαδραστικό χάρτη leaflet μέσω της βιβλιοθήκης folium
import folium
coords = list(wgs84_pt.coords)[0]
#Create the map
my_map = folium.Map(location = coords, zoom_start = 13)
# Add marker
folium.Marker(coords, popup = 'Volos').add_to(my_map)
#Display the map
my_map
Η βιβλιοθήκη Geopandas¶
Μέχρι τώρα είδαμε πως μπορούμε να διαχειριζόμαστε γεωμετρικά δεδομένα με την βιβλιοθήκη shapely.
Όμως τα γεωγραφικά δεδομένα δεν περιλαμβάνουν μόνο της γεωγραφική πληροφορία και την γεωμετρική τους δομή αλλά συνοδεύονται από μια σειρά περιγραφικών δεδομένων.
Η βιβλιοθήκη geopandas είναι μια βιβλιοθήκη της Python η οποία υποστηρίζει την ανάγνωση, επεξεργασία, ανάλυση και εγγραφή γεωγραφικών και παράλληλα περιγραφικών δεδομένων. Αποτελεί επέκταση της βιβλιοθήκης pandas και «κληρονομεί» τα χαρακτηριστικά και τις δυνατότητές της.
Στην υφιστάμενη δομή της pandas, η geopandas υποστηρίζει την γεωμετρία με την προσθήκη μιας νέας στήλης και το γεωγραφικό σύστημα αναφοράς.
Εισάγουμε την βιβλιοθήκη με τον παρακάτω τρόπο
import geopandas as gpd
Για να διαβάσουμε ένα αρχείο shapefile χρησιμοποιούμε την συνάρτηση gpd.read_file. Το αντικείμενο που προκύπτει είναι τύπου geopandas.geodataframe.GeoDataFrame
# Import shapefile using geopandas
dhmoi = gpd.read_file("../docs/dhmoi.gpkg")
type(dhmoi)
geopandas.geodataframe.GeoDataFrame
Η στήλη της γεωμετρίας ονομάζεται προκαθορισμένα geometry και είναι αντικείμενο geopandas.geoseries.GeoSeries που περιλαμβάνει αντικείμενα shapely.
type(dhmoi["geometry"])
geopandas.geoseries.GeoSeries
Με την μέθοδο geom_type μπορούμε να δούμε τον γεωμετρικό τύπο κάθε εγγραφής.
dhmoi.geom_type
0 MultiPolygon
1 MultiPolygon
2 MultiPolygon
3 MultiPolygon
4 MultiPolygon
...
320 MultiPolygon
321 MultiPolygon
322 MultiPolygon
323 MultiPolygon
324 MultiPolygon
Length: 325, dtype: object
Ας δούμε τις αρχικές εγγραφές του αρχείου
dhmoi.head()
| OBJECTID | X | Y | Name | CodeELSTAT | PopM01 | PopF01 | PopTot01 | UnemrM01 | UnemrF01 | UnemrT01 | PrSect01 | Foreig01 | Income01 | Perif | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 616259.007833 | 4.551127e+06 | KOMOTINIS | 0101 | 29967.0 | 31534.0 | 61501.0 | 4.726948 | 6.051873 | 5.221281 | 21.817852 | 1.996780 | 10969.315256 | 1 | MULTIPOLYGON (((631058.735 4574912.809, 631431... |
| 1 | 2 | 644783.135624 | 4.561364e+06 | ARRIANON | 0102 | 8952.0 | 9307.0 | 18259.0 | 1.955813 | 2.390181 | 2.154240 | 92.912436 | 0.038337 | 6398.641486 | 1 | MULTIPOLYGON (((661976.597 4554253.076, 662118... |
| 2 | 3 | 602100.950461 | 4.555689e+06 | IASMOU | 0103 | 7290.0 | 7561.0 | 14851.0 | 5.653884 | 6.671554 | 6.062390 | 70.206767 | 1.339977 | 6608.897469 | 1 | MULTIPOLYGON (((603623.489 4566376.502, 603657... |
| 3 | 4 | 633712.468442 | 4.539548e+06 | MARONEIAS - SAPON | 0104 | 8284.0 | 8342.0 | 16626.0 | 5.840957 | 10.270499 | 7.421934 | 57.748085 | 0.890172 | 6700.375345 | 1 | MULTIPOLYGON (((650998.525 4557469.860, 651446... |
| 4 | 5 | 517831.195188 | 4.572544e+06 | DRAMAS | 0201 | 28041.0 | 29326.0 | 57367.0 | 8.831689 | 14.176926 | 10.851115 | 7.298838 | 2.961633 | 11299.106498 | 2 | MULTIPOLYGON (((521986.700 4601855.354, 522271... |
dhmoi.plot()
<AxesSubplot:>
Χαρτογραφική απόδοση σε διαδραστικό χάρτη
dhmoi.explore(legend=False)